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ABSTRACT 

In this Letter, we investigate the claim of Valluri et al. (2003), namely that the use of 
Schwarzschild dynamical (orbit-based) models leads to an indeterminacy regarding the es- 
timation of the free parameters like the central black hole mass and the stellar mass-to-light 
ratio of the galaxy under study. We examine this issue with semi-analytic two-integral mod- 
els, which are not affected by the intrinsic degeneracy of three-integral systems. We however 
confirm the Valluri et al. result, and observe the so-called widening of the contours as 
the orbit library is expanded. We also show that, although two-dimensional data coverage 
help in constraining the orbital structure of the modelled galaxy, it does not in principle solve 
the indeterminacy issue, which mostly originates from the discretization of such an inverse 
problem. We show that adding regularization constraints stabilizes the confidence level con- 
tours on which the estimation of black hole mass and stellar mass-to-light ratio are based. 
We therefore propose to systematically use regularization as a tool to prevent the solution to 
depend on the orbit library. Regularization, however, introduces an unavoidable bias on the 
derived solutions. We hope that the present Letter will trigger some more research directed at 
a better understanding of the issues addressed here. 

Key words: galaxies: elliptical and lenticular galaxies: structure galaxies: nuclei galaxies: 
stellar dynamics 



1 INTRODUCTION 

A recent paper by Valluri, Merritt & Emsellem (2003; hereafter 
VME03) demonstrated that axisymmetric orbit-based modelling 
algorithms have only a limited ability to recover two fundamental 
parameters characterizing the total mass distribution, namely the 
stellar mass-to-light ratio T and the black hole (BH) mass M,. 
These authors showed that if the number of orbits is increased, the 
confidence contours defining the best fit values for [M. ; T] and 
their associated error bars widen up, hence significantly increasing 
the uncertainty of the BH measurement (see e.g., their Figure 3). 
In the case of M32, they showed that a range of models with M, 
values between 1 x 10® Mq and 6 x 10® Mq provide equally good 
fit to the data discussed by van der Marel et al. (1998). This con- 
trasts with the claim of other groups, e.g., "the range of black hole 
mass uncertainties is from 5% to 70%, with an average uncertainty 
around 20 %" (Gebhardt 2003). 

We have investigated this question using component-based 
f[E, L,) models, where f[E, Lz) is the distribution function (DF) 
depending only on the energy E and the vertical angular momen- 
tum Lz (Sect.|2j. We confirm the same widening effect reported by 
VME03, albeit with two-integral models instead of three-integral 



models. In Sect. 12. 31 we argue that the specific effect we describe 
here is due to the discretization of phase-space (i.e., the sampling 
into a finite number of orbits or components) and is therefore not 
related to the peculiar form of the distribution function. 

We illustrate our claim within a realistic context, namely using 
the luminous mass model, as well as the SAURON and OASIS ^ 
setups (Bacon et al. 2001, de Zeeuw et al. 2002) for NGC 3377 (see 
Copin et al. 2003, hereafter C03): these are described in Sect.|2| 

We have considered two options in order to examine the relia- 
bility of the T and M, determinations. First, we have investigated 
whether the spatial extension of the data helps to better constrain 
the model (Sect.|3). Second, we have included regularization in the 
fitting process (Sect.|4j and show that it provides a way to stabilize 
the widening effect in the two-integral case. 
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SAURON and OASIS are integral-field spectrographs 
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Figure 1. x'^-contours without any regularization. Models including only the OASIS constraints are shown on the top row, whereas the bottom row displays 
models with only the SAURON constraints. The number of components increases from left to right (see insets). The third, thicker, contour (counting from the 
innermost one) corresponds to the formal 99.7% confidence level with two degrees of freedom. 



2 INPUT CONSTRAINTS AND DYNAMICAL 
MODELLING 

2.1 A two-integral model 

One of the advantages of working with a two-integral DF is unique- 
ness: there is a unique f{E, Lz) DF that matches simultaneously 
the mass density, which is related to the even part of the DF, and 
the mean velocity field which is fully specified by the odd part of 
the DF (see e.g., Merritt 1996b). There is therefore no intrinsic de- 
generacy in recovering the input parameters [A/.; T] for the two- 
integral models. Note, however, that in the three-integral case the 
model has some additional freedom to adjust its orbital anisotropy 
(e.g., in the three-integral case, the radial velocity dispersion can 
be different than the vertical one). Thus, three-integral models with 
[M, \ T] values very different from the true input values may still 
be found to fit the data equally well (see VME03). Because of this 
extra freedom, the confidence level contours will always be larger 
in the three-integral case. 

Another advantage of the two-integral models is that we can 
easily check directly their entire dynamical structure by comparing 
the reconstructed DF with the true input DF. In our view, these two 
properties make the two-integral case a genuine test case. 



2.2 The mass model and kinematic constraints 

For the luminous density distribution we use our Multi-Gaussian 
Expansion (MGE) model of NGC 3377 (see C03, and Emsellem 
et al. 1994). The stellar T/ has been fixed to 2.4 and a dark com- 
ponent in the form of a 10* Mq black hole has been included at 
the center: this fixes the gravitational potential of our axisymmetric 
model. 

Except for a few notable cases (see e.g., Cappellari et al. 2002; 
Verolme et al. 2002, C03), the vast majority of BH determinations 
are based on one-dimensional long-slit spectra, leaving most of the 
(x' , y') sky plane unconstrained. One might reasonably expect that 
data densely covering the projected plane of the sky (x' , y') would 
help to constrain the models. This is indeed emphasized in Verolme 
et al. (2002) and in C03, where it is shown that the estimates of the 
free parameters of the models (e.g., inclination, M,, T) are signif- 
icantly better constrained when two-dimensional photometric and 
kinematic data are included. We have therefore tested two different 
regimes, namely with two-dimensional constraints having different 
spatial extensions (see C03; Bacon et al. 2001): 

(i) a central dataset, with the characteristics of the OASIS setup 
covering the inner 4" by 2" and. 
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Figure 2. Same as figurelTlbut for models including regularization 



(ii) a spatially extended dataset, with the full SAURON field of 
view covering 46" by 42" . 

From the MGE mass model, we have generated artificial kine- 
matic constraints assuming a semi-isotropic distribution function 
f{E,L^) using the Hunter-Qian (1993) method (HQ hereafter). 
The even part of such a distribution function (DF) is uniquely spec- 
ified by the MGE mass model, and we chose the odd part so as 
to roughly mimic the observed SAURON and OAS I S kinematics of 
NGC 3377 (C03). The observables of the HQ model (i.e., the line- 
of-sight velocity profiles, hereafter VPs) are then computed for 
both setups: they are convolved with the appropriate point spread 
functions and binned into the respective spatial elements, and will 
be used to constrain the dynamical models^ . 

2.3 The f(E, Lz ) component-based Schwarzschild models 

We now turn to the building of a Schwarzschild dynamical model 
constrained by the fake photometric and kinematic datasets de- 
scribed in Sect. 12.21 We use the f{E, L^) component-based mod- 
els of Cretton et al. 1999 (hereafter C99). These authors described 

^ In the remainder of this paper, "SAURON constraints" refer to the ailificial 
data derived from the HQ method, using the SAURON setup as in C03 (pixel 
size. Point Spread Function, etc). 



how to derive two-integral models (with or without central black 
holes) using only these f{E, Lz) components. These models were 
fully specified by reconstruction of the corresponding DFs (see 
their Figs. 6, 7 and 8), which were compared to the true under- 
lying DFs obtained via the HQ method (see Verolme & de Zeeuw 
2002 for similar computations). Following C99, and using the grav- 
itational potential corresponding to the MGE mass model (Sect.|2j, 
we have created a library of two-integral components by sampling 
the energy E with N e values and, at each energy, by sampling the 
angular momentum with Nl^ values. The final solutions are ob- 
tained with a non-negative least square algorithm (NNLS, see C99 
and Lawson & Hanson 1974 for details). 

In the next Sections, we show that the widening effect reported 
by VME03 is also present with two-integral models, despite the 
absence of an intrinsic degeneracy for such models. 



3 SPATIAL EXTENSION 

FigureQshows the contours obtained without regularization us- 
ing 1. only the OASIS constraints, covering the central parts (top 
row), or 2. only the SAURON constraints with a large field of view 
(bottom row). 

The OAS I S contours exhibit a clear correlation between the 
M, and T / : a higher black hole mass needs to be compensated by 
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Figure 3. Distributions functions: tlie dotted line is the "true" 
DFiiq{E , Lz) derived using tlie HQ algoritlim and tlie full line is the re- 
constructed DFnnls(-E, Lz) of the model with [A/.; T/]=[10'*; 2.4]. 

a lower T/. With their spatial extension, the SAURON-like observ- 
ables provide a much better constraint on T/, but obviously fail to 
restrict the allowed range for M, . This is an expected behaviour for 
such models, as discussed in details in C03. What is more surpris- 
ing is that for both cases, when the number of components Ncomp 
is increased, a larger region of the [M,; T i] plane is enclosed by 
the 99.7% confidence level (Scr error in two degrees of freedom, 
which corresponds to Ax^ = 11.8, see Press et al. 1992). We have 
increased Ncomp by up to a factor of 16 and still, we see no sign of 
convergence. Using our two-integral models, we therefore see an 
effect similar to the one reported by VME03. 

The wider field of view does not seem to stop the x^-contours 
from getting broader; as emphasized by the bottom row of Fig. □ 
it merely slows the effect down. However for the SAURON setup, 
we can still observe a significant increase in the allowed range of 
T I (i.e., within the 3a contour), particularly when N_e = 80 and 
Ni, = 28. 



4 TOWARDS REGULARIZATION 

In the absence of regularization, the weights of many two-integral 
components are set to zero by the NNLS algorithm: only Nconst 
components can have a non zero weight, where Nconst is the num- 
ber of constraints (see C99). Such a behaviour has been emphasized 
by e.g., Cappellari et al. (2002), and Verolme & de Zeeuw (2002). 

By increasing Ncomp one improves the sampling of the inte- 
gral space, which also corresponds to a better spatial sampling. In 
this sense, with a bigger library NNLS has more flexibility to fit the 
(same) data. An equally good fit to the data can then be obtained 
using a different model than the true input one. A larger fraction of 
the models in [Af, ; T]-space are now acceptable, because NNLS 
has more freedom to select the appropriate orbits. The discretiza- 
tion of the phase space implies that the problem of recovering the 



orbital weights become ill-conditioned. We believe the same effect 
occurs with three-integral orbit-based models and suggest that it is 
partly at the root of the contour widening problem (VME03). 

Including regularization forces all the components to have a 
non zero weight and therefore greatly reduces (or eliminates) this 
freedom. Furthermore, with such jagged weight distributions, the 
solutions without regularization are definitely far away from phys- 
ical solutions. We do not claim that our way of performing regu- 
larization is optimal or represent in any sense what is taking place 
in real galaxies. As emphasized in Sect. |3|and in VEM03, unreg- 
ularized solutions will depend on the size of the orbit library used 
in the analysis. However, regularization introduces an unavoidable 
bias in the solutions the algorithm will find. 

Regularization is done here in the same way as described in 
C99, i.e., by minimizing the second derivatives of the weights on a 
grid in Integral space. 



4.1 Results 

Figures |2| shows the results (as contours) of using the OASIS 
(top row) or SAURON (bottom row) constraints, when regulariza- 
tion is included. The amount of regularization was set such as to 
optimally reproduce the input HQ-DF (see Fig. 8 of C99). Some 
other schemes to tune the optimal regularization parameter exist 
(e.g., generalized cross-validation, Wahba 1990; Merritt 1997). 

In both SAURON and OASIS setups, the contours are stable in 
the sense that, while increasing Ncomp by up to a factor 16, they 
do not show the widening effect observed without regularization. 
This can be intuitively understood by examining the size of the 
matrix problem which NNLS is trying to solve. Without regular- 
ization, the input matrix containing the observables for each com- 
ponent has Ncomp columns (number of components) and Nconst 
lines (number of constraints). When Ncomp is increased reach- 
ing values significantly larger than Nconst, the problem becomes 
strongly underconstrained, and most of the resulting weights are 
set to zero (see above). When regularization is introduced, the num- 
ber of columns stays constant, while the number of lines becomes 
Ncomp + Nconst X Nint ( Nint is the uumbcr of integral of motions, 
equal to 2 in the present two-integral case). As Ncomp is increased, 
the numbers of columns and lines of the matrix increase simultane- 
ously, their ratio converging to 1 /Nmt . Although it depends on the 
regularization scheme used, these additional constraints can thus 
qualitatively change the nature of the inverse problem. 

We should emphasize here that the results presented in this 
Section are no proof that the contours will not get wider for a much 
larger Ncomp. Nevertheless, it is clearly an improvement over the 
unregularized case, the results of which depend on the size of the 
orbit library. 



4.2 Distribution functions 

A more conclusive test is to examine the DF corresponding to the 
weights recovered by the least-square routine. Figure |3| compares 
the reconstructed DFnnls obtained from the NNLS fit to the input 
DFhq derived from HQ. For each value of the Energy, Lz runs 
monotonically from I/z,min to Lz,max, covering Nl^ (=7 here) 
components. The next component corresponds to the next value 
of the energy, etc. Overall, the agreement is good, except for dis- 
crepancies at the peaks near Z/z,max- A similar test presented by 
C99 on a fake model of M 32 showed a significantly better fit (with 
RMS residuals of the order of a few percents). The difference lies 
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certainly in the contribution of a disk-like component in the fake 
model of NGC 3377, which would probably require a better sam- 
pling in Lz- It is likely that some more advanced regularization 
scheme could improve the recovery of these steep gradients in the 
DF, but this is outside the scope of the present paper. 



5 COMPARISON WITH OTHER MODELS 

Gebhardt (2003) claimed that increasing the number of orbits does 
not produce the contour widening effect described in VME03 and 
in the present paper. The cause for this discrepancy is, at the mo- 
ment, unclear to us. 

Models with and without regularization look surprisingly si- 
milar when N_B=20 and Nl,=7, which are typical values used in 
many orbit-based dynamical models (see e.g., vdM98, C99, Cret- 
ton & van den Bosch 1999, Verolme et al. 2002, Cappellari et 
al. 2002). We suggest that this coincidence has (at least) two con- 
sequences. Firstly, this could trigger the misconception that reg- 
ularization does not affect the indeterminacy in the estimation of 
[M,; T], since the contours (with or without regularization) are 
similar. Secondly, this could explain the good apparent agreement 
sometimes found between regularized and unregularized codes. If 
this comparison had been performed using a much larger number 
of orbits (in the unregularized case), we predict that a different con- 
clusion would have been reached. 



6 CONCLUSION 

In this Letter, we have used two-integral models to test the claim 
made by VME03 that contours derived using the Schwarzschild 
orbit-based technique widen as the orbit library is expanded. 

• We confirm this effect with two-integral models for which the 
intrinsic degeneracy discussed in VME03 does not exist. 

• We show that the widening of the contours occurs even 
when we use two-dimensional kinematics constraints. 

• We suggest that the origin of such an effect partly lies in the 
discretization of phase space. 

• We finally propose to systematically use regularization as a 
way to stabilize the diagrams produced by such techniques, so 
that the obtained solution does not significantly depend on the input 
orbit library. 

Regularization schemes have been often advocated to derive 
stable solutions of ill-conditioned inverse problems (e.g., Golub et 
al. 1999, and references therein; Merritt 1996b). The main worry 
associated with regularization is that it biases the solution found 
by the fitting algorithm (see e.g., Merritt, 1996a and Dehnen, 2001 
for similar questions in N-body codes). This bias may indeed not 
correspond to the true characteristics of the object under study. The 
input two-integral DF we chose for the present test is, by construc- 
tion, a smooth (semi-analytical) function of the energy and angular 
momentum. Our analysis shows that in this case a simple regular- 
ization scheme with the right choice of smoothing parameter can be 
applied to recover the full (two-integral) DF reasonably well (see 
also C99 for a similar demonstration on a less complex DF). When 
modelling real galaxies, however, there is no obvious way to select 
the proper level of smoothing, an incorrect choice leading to either 
degenerate solutions (if undersmoothed; see Fig.0 or biased solu- 
tions (if oversmoothed). The sampling of the orbit library may be 



further optimized such as to ensure a relatively homogeneous cov- 
erage of the presumed components of the galaxy (core, disk, bulge, 
etc). As mentioned in Sect. 14.21 it is likely that an adaptive regu- 
larization scheme can be designed, e.g., taking into account some a 
priori knowledge (or preconception) we have on a galaxy (e.g., the 
presence of a bright disk). 

The stability of regularized solutions however comes at the 
cost of an unavoidable bias. The resulting best fit parameters and 
confidence range for [A/,; T] in a galaxy may not necessarily re- 
flect their true values. We emphasize that, in this Letter, we only 
treat the indeterminacy caused by the discrete nature of the in- 
verse problem and not the additional degeneracy associated with 
3-integral models, as reported by VME03. We therefore hope that 
the present Letter will trigger further studies to address these issues. 
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